###### Title ##########################
## Insecurity and Support for Female Leadership in Conflict States: 
## Evidence from Afghanistan by Jasmine Bhatia and
## Steve L. Monroe.
## Replication Command File  
## August 2023
##

# This file contains analysis and tables of the conjoint survey experiment
# presented in both the main article and supplementary information


### Table of Contents  ##################################
## 1. Load packages
## 2. Upload data, clean variables
## 3. Balance Tables
## 4. H1: Insecurity and Support for Female Leadership
## 5. H2: Insecurity and Support for Female Leadership Across Genders
## 6. Mechanisms


## Note: All code run using both R version 3.6.3 (2020-02-29)

######## 1. Load Necessary Packages ##############################

library(cjoint)
library(dplyr)
library(tidyr)
library(xtable)
library(tidyverse)
library(stats)
library(survey)
library(Rcpp)
library(estimatr)
library(knitr)
library(DT)
library("here")
library(cli)
library(devtools)
devtools::install_github("m-freitag/cjpowR")
library(cjpowR)
library(stargazer)
library(MASS)
library(mlogit)
library("nnet")
library("scales")

# set working directory and load dataset

getwd()

### 2. Load Data and Clean Variables #######################
##     


data <- readRDS("2019_Afghan_data_R.RDS")


## Select relevant variables for analysis ###
## These are variables pertaining to insecurity (IV), gender equality (DV 1),
## support for the military (DV 2), and control variables.


data2 <- data[, c("x15", "x15b", "x16", "x17_9", "x17_11", 
                  "x34a", "x34b", "x35a", "x35d",
                  "g8a", "x26g", "x35b",
                  "g8c", "x35f", "x35e",
                  "x69", "x378", "x440",
                  "z1", "z2", "z6", "z9", "z10", "z13b",
                  "edu", "m4", "wave", "m7", "m1",
                  "g6b", "x280", "x37bg", "x41d",
                  "x37b", "x70a", "x70b" 
                  )]


colnames(data2) <- c("fear_frequency_a", "fear_frequency_b", "violence_experience", "violence_insurgents", "violence_army",
                     "confidence_ana", "confidence_police", "fair_ana", "security_improve_ana",
                     "ana_security_impact", "fear_ana", "ana_unprofessional",
                     "security_impact_police", "protects_civilians_ana", "ana_professional",
                     "leadership_gender", "Taliban_negotiations_good" , "Taliban_peace_support",
                     "gender", "age", "edu_self_reported", "married",
                     "ethnicity", "women_contribute_income", 
                     "edu", "region", "wave", "province", "respondent_id",
                     "rate_security_area_1", "rate_security_area_2", "gov_handling_security", "women_indv_right",
                     "prov_gov", "women_parl", "women_prov")


## Clean variables 

# Independent Variable

data2$female <- ifelse(data2$gender == "Female", 1, 0)

### Insecurity Variables 

## How often do you feel How often do you fear for your own personal safety 
#or security or for that of your family these days? Often, sometimes, rarely, or never?

data2$fear_numeric_a  <- case_when(
  data2$fear_frequency_a == "Often" ~ 4,
  data2$fear_frequency_a == "Sometimes" ~ 3,
  data2$fear_frequency_a == "Rarely" ~ 2,
  data2$fear_frequency_a == "Never" ~ 1,
  TRUE ~ NA_real_
)

data2$fear_numeric_no_na_a <- data2$fear_numeric_a

data2$fear_numeric_no_na_a[is.na(data2$fear_numeric_no_na_a)] <- 0


data2$fear_numeric_b  <- case_when(
  data2$fear_frequency_b == "Always" ~ 4,
  data2$fear_frequency_b == "Often" ~ 4,
  data2$fear_frequency_b == "Sometimes" ~ 3,
  data2$fear_frequency_b == "Rarely" ~ 2,
  data2$fear_frequency_b == "Never" ~ 1,
  TRUE ~ NA_real_
)

data2$fear_numeric_no_na_b <- data2$fear_numeric_b

data2$fear_numeric_no_na_b[is.na(data2$fear_numeric_no_na_b)] <- 0


### Fear Dummy that incorporates both

data2$fear_numeric <- data2$fear_numeric_no_na_a + data2$fear_numeric_no_na_b

### Remove 0 observations (these would be Refused to answers = 541)

data2$fear_numeric <- ifelse(data2$fear_numeric > 0, data2$fear_numeric,
                             NA)


# Equals one if feel fear sometimes or always

data2$fear_dummy <- ifelse(data2$fear_numeric > 2, 1, 0)


# One if feels fear always
data2$fear_dummy_high <- ifelse(data2$fear_numeric > 3, 1, 0)

# Have you or has anyone in your family been a victim of violence 
# or of some criminal act in your home or community in the past year? (1 = Yes)


data2$violence_experience_dummy  <- case_when(
  data2$violence_experience == "yes" ~ 1,
  data2$violence_experience == "no" ~ 0,
  TRUE ~ NA_real_
)



## Dependent Variable: 1 Support for Female Political Representation  

## Should be political leadership positions be 1) Mostly for women, 2) Equal, 3) Mostly for Men

data2$female_support_continuous <- ifelse(data2$leadership_gender == "mostly for women", 1,
                       ifelse(data2$leadership_gender == "equal for both men and women", 0,
                        ifelse(data2$leadership_gender == "anyone based on merit", 0,      
                        ifelse(data2$leadership_gender == "mostly for men", -1,
                        ifelse(data2$leadership_gender == "women should do house work", -1,      
                        NA)))))


data2$female_support_unordered <- ifelse(data2$leadership_gender == "mostly for women", "mostly for women",
                                  ifelse(data2$leadership_gender == "equal for both men and women", "equal",     
                                  ifelse(data2$leadership_gender == "anyone based on merit", "equal",
                                   ifelse(data2$leadership_gender == "mostly for men", "men",
                                          NA))))      
                                         
                                         
                                  

data2$female_support_ordered <- factor(data2$female_support_continuous,
                                       levels= c("-1", "0", "1"))

data2$female_support_continuous2 <- ifelse(data2$leadership_gender == "mostly for women", 1,
                                     ifelse(data2$leadership_gender == "equal for both men and women", 0,
                                     ifelse(data2$leadership_gender == "anyone based on merit", NA,      
                                       ifelse(data2$leadership_gender == "mostly for men", -1,
                                      ifelse(data2$leadership_gender == "women should do house work", -1,      
                                        NA)))))

data2$mostly_men <- ifelse(data2$leadership_gender == "mostly for men", 1, 
                           ifelse(data2$leadership_gender == "women should do house work", 1,      
                                  0))

data2$mostly_men2 <- ifelse(data2$leadership_gender == "mostly for men", 1, 
                           ifelse(data2$leadership_gender == "women should do house work", 1,
                          ifelse(data2$leadership_gender == "anyone based on merit", NA,
                          0)))      


data2$gender_equality <- ifelse(data2$leadership_gender == "equal for both men and women", 1, 0)

data2$gender_equality2 <- ifelse(data2$leadership_gender == "equal for both men and women", 1,
                          ifelse(data2$leadership_gender == "anyone based on merit", NA,      
                                 0))


data2$mostly_women <- ifelse(data2$leadership_gender == "mostly for women", 1, 0)

data2$mostly_women2 <- ifelse(data2$leadership_gender == "mostly for women", 1, 
                      ifelse(data2$leadership_gender == "anyone based on merit", NA,   
                         0))

## Examples of Patriarchy

# Are you opposed to women representing you in parliament

data2$women_parl_binary <- ifelse(data2$women_parl == "Yes", 1,
                           ifelse(data2$women_parl == "No", 0,
                           NA))

# Are you opposed to women representing you in provincial gov?
data2$women_prov_binary <- ifelse(data2$women_prov == "Yes", 1,
                            ifelse(data2$women_prov == "No", 0,
                             NA))


## Mechanism: Support for prov gov

data2$prov_gov_support_rank <- ifelse(data2$prov_gov == "very bad job", 0,
                               ifelse(data2$prov_gov == "somewhat bad job", 1,        
                               ifelse(data2$prov_gov == "somewhat good job", 2,       
                              ifelse(data2$prov_gov == "very good job", 3, 
                              NA))))

data2$prov_gov_support_dummy  <- ifelse(data2$prov_gov_support_rank > 1, 1, 0)    



## Control Variables

data2$married_dummy <- ifelse(data2$married == "married", 1, 0)

data2$wave_factor <- as.factor(data2$wave)

data2$edu_numeric <- ifelse(data$edu == "No formal school", 0,
                     ifelse(data$edu == "Informal schooling at home or literacy class", 0,
                     ifelse(data$edu == "Primary school (classes 1-6, complete or incomplete)", 1,
                     ifelse(data$edu == "Secondary school (classes 7-9, complete or incomplete)", 2,
                     ifelse(data$edu == "High school (classes 10-12, complete or incomplete)", 3,
                     ifelse(data$edu == "University education (classes 13-18, complete or incomplete)", 4,
                     NA))))))




# Variable for whether respondents lived in a province where 
# the conjoint survey experiment was carried out

data2$experiment <- ifelse(data2$province == "Kunduz", 1,
                    ifelse(data2$province == "Balkh", 1,      
                    ifelse(data2$province == "Sar-e-pul", 1,
                      0)))       


# Dataset for respondents in surveys where we did not run survey experiment

data3b <- data2[data2$experiment == 0,]


## Summary Statistics about men vs. women and their leadership preferences


## Do men have low levels of support for women? (floor)

male <- data2[data2$female == 0, ]

female <- data2[data2$female == 1, ]


## 3. Descriptive Statistics #####

## SI Section 3.2

# SI Table 18 (hand code)


summary(female[female$violence_experience_dummy == 0,]$mostly_men2)

sd(female[female$violence_experience_dummy == 0,]$mostly_men2, na.rm = TRUE)

summary(female[female$violence_experience_dummy == 1,]$mostly_men2)

sd(female[female$violence_experience_dummy == 1,]$mostly_men2, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 0,]$mostly_men2)

sd(male[male$violence_experience_dummy == 0,]$mostly_men2, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 1,]$mostly_men2)

sd(male[male$violence_experience_dummy == 1,]$mostly_men2, na.rm = TRUE)


## SI Table 19 (hand code)


summary(female[female$violence_experience_dummy == 0,]$mostly_women2)

sd(female[female$violence_experience_dummy == 0,]$mostly_women2, na.rm = TRUE)

summary(female[female$violence_experience_dummy == 1,]$mostly_women2)

sd(female[female$violence_experience_dummy == 1,]$mostly_women2, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 0,]$mostly_women2)

sd(male[male$violence_experience_dummy == 0,]$mostly_women2, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 1,]$mostly_women2)

sd(male[male$violence_experience_dummy == 1,]$mostly_women2, na.rm = TRUE)


# SI Table 20 Pro-Women Score

summary(female[female$violence_experience_dummy == 0,]$female_support_continuous)

sd(female[female$violence_experience_dummy == 0,]$female_support_continuous, na.rm = TRUE)

summary(female[female$violence_experience_dummy == 1,]$female_support_continuous)

sd(female[female$violence_experience_dummy == 1,]$female_support_continuous, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 0,]$female_support_continuous)

sd(male[male$violence_experience_dummy == 0,]$female_support_continuous, na.rm = TRUE)

summary(male[male$violence_experience_dummy == 1,]$female_support_continuous)

sd(male[male$violence_experience_dummy == 1,]$female_support_continuous, na.rm = TRUE)


### Figure 4 in Manuscript

## Subset secure respondents

secure <- data2[data2$violence_experience_dummy == 0, ]


secure_mean <- as.vector(c(mean(secure$mostly_men2, na.rm = TRUE), 
                             mean(secure$gender_equality2, na.rm =TRUE),
                             mean(secure$mostly_women2, na.rm = TRUE)))



## Subset insecure respondents
insecure <- data2[data2$violence_experience_dummy == 1, ]


insecure_mean <- as.vector(c(mean(insecure$mostly_men2, na.rm = TRUE), 
            mean(insecure$gender_equality2, na.rm =TRUE),
            mean(insecure$mostly_women2, na.rm = TRUE)))
            
    
table <- c(secure_mean, insecure_mean)

answer <- c("Mostly Men", "Equal", "Mostly Women", "Mostly Men", "Equal", "Mostly Women")

group <- c("Secure", "Secure", "Secure", 
           "Insecure", "Insecure", "Insecure")

table <- cbind(answer, table, group)

graph <- as.data.frame(table)

graph$table <- as.numeric(graph$table)

graph$answer <- factor(graph$answer, order = TRUE, 
                       levels =c("Mostly Men", "Equal", "Mostly Women"))

graph$group <- factor(graph$group, order = TRUE,
                      levels = c("Secure", "Insecure"))

figure4 <-ggplot(graph, aes(fill = group, x=answer, y=table)) +
              geom_bar(position= "dodge", 
              stat="identity", width=0.5, color = "black")+
              scale_y_continuous(labels = number_format(accuracy = 0.1))+
              scale_fill_manual(values = c("#CCCCCC", "#000000"))+
              theme(
              axis.title.x = element_blank(),
              axis.title.y = element_blank(),
              legend.title=element_blank())
              

##### 4. Balance Tables between Insecure and Secure Respondents ####


## SI Table 15

balance_table_asia <- data2 %>%
  group_by(violence_experience_dummy) %>%
  summarize(
    avg_female = mean(female, na.rm=TRUE),
    sd_female = sd(female, na.rm=TRUE),
    avg_age = mean(age, na.rm=TRUE),
    sd_age = sd(age, na.rm=TRUE),
    avg_edu = mean(edu_numeric, na.rm=TRUE),
    sd_edu = sd(edu_numeric, na.rm=TRUE),
    avg_marriage = mean(married_dummy, na.rm=TRUE),
    sd_marriage = sd(married_dummy, na.rm=TRUE),
    avg_experiment = mean(experiment, na.rm = TRUE),
    sd_experiment = sd(experiment, na.rm = TRUE),
    avg_prov_support = mean(prov_gov_support_dummy, na.rm = TRUE),
    sd_prov_support = sd(prov_gov_support_dummy, na.rm = TRUE),
    count = n())


pvl <- data2 %>%
  summarise(
    pval_ = t.test(female~ violence_experience_dummy, 
                   var.equal = TRUE)$p.value,
    pval_age = t.test(age ~ violence_experience_dummy,
                      var.equal = TRUE)$p.value,
    pval_edu = t.test(edu_numeric ~ violence_experience_dummy, 
                      var.equal = TRUE)$p.value,
    pval_marriage = t.test(married_dummy ~ violence_experience_dummy, 
                           var.equal = TRUE)$p.value,
    pval_experiment = t.test(experiment ~ violence_experience_dummy, 
                             var.equal = TRUE)$p.value,
    pval_prov_support = t.test(prov_gov_support_dummy ~ violence_experience_dummy, 
                               var.equal = TRUE)$p.value,
  )




control <- as.vector(t(balance_table_asia[1,-c(1,14)]))
treatment <- as.vector(t(balance_table_asia[2,-c(1,14)]))


p_val_one_asia <- as.vector(unlist(c(pvl[1], NA, pvl[2], NA, pvl[3], NA, 
                                     pvl[4], NA, pvl[5], NA, pvl[6], NA
)))


b_vars <- c(
  "Female", "",
  "Age", "",
  "Education", "",
  "Marital Status", "",
  "Survey Experiment Province", "",
  "Prov Support", "",
  "Respondents")


counts <- balance_table_asia$count

cb_tab <- data.frame(matrix(nrow = 13, ncol = 4))

cb_tab[,1] <- c(b_vars)
cb_tab[,2] <- c(control, counts[1])
cb_tab[,3] <- c(treatment, counts[2])
cb_tab[,4] <- c(p_val_one_asia, NA)


# Round all numbers to 3 digits
cb_tab <- cb_tab %>%
  purrr::modify_if(~is.numeric(.), ~round(., 3))


# Add descriptive row at top
cb_tab1_asia<- data.frame(rbind(c("", "Mean \n (SD)", "Mean \n (SD)", 
                                  "P Value"),
                                cb_tab))

# Add column names
names(cb_tab1_asia) <- c("", "No Violence", "Violence", "")

View(cb_tab1_asia) # SI Table 15 Balance 


## SI Table 16: Exposure to Insecurity (Female Respondents)


balance_table_women <- female %>%
  group_by(violence_experience_dummy) %>%
  summarize(
    avg_age = mean(age, na.rm=TRUE),
    sd_age = sd(age, na.rm=TRUE),
    avg_edu = mean(edu_numeric, na.rm=TRUE),
    sd_edu = sd(edu_numeric, na.rm=TRUE),
    avg_marriage = mean(married_dummy, na.rm=TRUE),
    sd_marriage = sd(married_dummy, na.rm=TRUE),
    avg_prov_support = mean(prov_gov_support_dummy, na.rm = TRUE),
    sd_prov_support = sd(prov_gov_support_dummy, na.rm = TRUE),
    count = n())


pvl <- female %>%
  summarise(
    pval_age = t.test(age ~ violence_experience_dummy,
                      var.equal = TRUE)$p.value,
    pval_edu = t.test(edu_numeric ~ violence_experience_dummy, 
                      var.equal = TRUE)$p.value,
    pval_marriage = t.test(married_dummy ~ violence_experience_dummy, 
                           var.equal = TRUE)$p.value,
    pval_prov_support = t.test(prov_gov_support_dummy ~ violence_experience_dummy, 
                               var.equal = TRUE)$p.value,
    
  )

control <- as.vector(t(balance_table_women[1,-c(1,10)]))
treatment <- as.vector(t(balance_table_women[2,-c(1,10)]))


p_val_one_asia <- as.vector(unlist(c(pvl[1], NA, pvl[2], NA, pvl[3], NA,
                                     pvl[4], NA
)))


b_vars <- c(
  "Age", "",
  "Education", "",
  "Marital Status", "",
  "Provincial Gov Support", "",
  "Respondents")


counts <- balance_table_women$count

cb_tab <- data.frame(matrix(nrow = 9, ncol = 4))

cb_tab[,1] <- c(b_vars)
cb_tab[,2] <- c(control, counts[1])
cb_tab[,3] <- c(treatment, counts[2])
cb_tab[,4] <- c(p_val_one_asia, NA)


# Round all numbers to 3 digits
cb_tab <- cb_tab %>%
  purrr::modify_if(~is.numeric(.), ~round(., 3))


# Add descriptive row at top
cb_tab2_asia<- data.frame(rbind(c("", "Mean \n (SD)", "Mean \n (SD)", 
                                  "P Value"),
                                cb_tab))

# Add column names
names(cb_tab2_asia) <- c("", "No Violence", "Violence", "")

View(cb_tab2_asia) # Balance Table  


## SI Table 17: Balance Table Exposure to Insecurity (Men) 


balance_table_men <- male %>%
  group_by(violence_experience_dummy) %>%
  summarize(
    avg_age = mean(age, na.rm=TRUE),
    sd_age = sd(age, na.rm=TRUE),
    avg_edu = mean(edu_numeric, na.rm=TRUE),
    sd_edu = sd(edu_numeric, na.rm=TRUE),
    avg_marriage = mean(married_dummy, na.rm=TRUE),
    sd_marriage = sd(married_dummy, na.rm=TRUE),
    avg_prov_support = mean(prov_gov_support_dummy, na.rm = TRUE),
    sd_prov_support = sd(prov_gov_support_dummy, na.rm = TRUE),
    count = n())


pvl <- male %>%
  summarise(
    pval_age = t.test(age ~ violence_experience_dummy,
                      var.equal = TRUE)$p.value,
    pval_edu = t.test(edu_numeric ~ violence_experience_dummy, 
                      var.equal = TRUE)$p.value,
    pval_marriage = t.test(married_dummy ~ violence_experience_dummy, 
                           var.equal = TRUE)$p.value,
    pval_prov_support = t.test(prov_gov_support_dummy ~ violence_experience_dummy, 
                               var.equal = TRUE)$p.value,
  )

control <- as.vector(t(balance_table_men[1,-c(1,10)]))
treatment <- as.vector(t(balance_table_men[2,-c(1,10)]))


p_val_one_asia <- as.vector(unlist(c(pvl[1], NA, pvl[2], NA, pvl[3], NA,
                                     pvl[4], NA
)))


b_vars <- c(
  "Age", "",
  "Education", "",
  "Marital Status", "",
  "Provincial Gov Support", "",
  "Respondents")


counts <- balance_table_men$count

cb_tab <- data.frame(matrix(nrow = 9, ncol = 4))

cb_tab[,1] <- c(b_vars)
cb_tab[,2] <- c(control, counts[1])
cb_tab[,3] <- c(treatment, counts[2])
cb_tab[,4] <- c(p_val_one_asia, NA)


# Round all numbers to 3 digits
cb_tab <- cb_tab %>%
  purrr::modify_if(~is.numeric(.), ~round(., 3))


# Add descriptive row at top
cb_tab3_asia<- data.frame(rbind(c("", "Mean \n (SD)", "Mean \n (SD)", 
                                  "P Value"),
                                cb_tab))

# Add column names
names(cb_tab3_asia) <- c("", "No Violence", "Violence", "")

View(cb_tab3_asia) # SI Table 17 



##### 5. H1: Insecurity and Support for Female Leadership #####

## SI Table 21 (Section 3.3)


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = data2) 


reg.1x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               female + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = data2, Hess=TRUE)


reg.2x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric +
                   married_dummy + province + wave_factor, 
                    data=data2)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3) # to get observations


# SI Table 22: Alternative Measure of Insecurity: Fear Dummy  


reg.1 <- lm(female_support_continuous ~ fear_dummy + female + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = data2) 


reg.1x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1



# See with ordinal

reg.2 <- polr(female_support_ordered  ~ fear_dummy + 
               female + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = data2, Hess=TRUE)


reg.2x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ fear_dummy +
                   female + age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=data2)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


# SI Table 23: No Wave fixed effects


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province,
            data = data2) 


reg.1x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1



# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               female + age + edu_numeric +
               married_dummy + province,
             data = data2, Hess=TRUE)


reg.2x <- na.omit(data2[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric +
                   married_dummy + province, 
                 data=data2)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



## SI Table 24:  Only look at three provinces where we ran
## survey experiment (Kunduz, Sar-e-pul, Balkh)

# kunduz, Sar-e-pul, balkh

kunduz <- data2[data2$province == "Kunduz",]
balkh <- data2[data2$province == "Balkh",]
sar <- data2[data2$province == "Sar-e-pul",]

data3 <- rbind(balkh, kunduz, sar)


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = data3) 


reg.1x <- na.omit(data3[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               female + age + edu_numeric +
               married_dummy + wave_factor,
             data = data3, Hess=TRUE)


reg.2x <- na.omit(data3[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric + 
                   married_dummy + province + wave_factor, 
                 data=data3)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



### SI Table 25: Insecurity and Support for Female Leadership (Hawks) 

hawks <- data2[data2$Taliban_peace_support != "Strongly support", ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province,
            data = hawks) 


reg.1x <- na.omit(hawks[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1




# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
                female + age + edu_numeric +
                married_dummy + province,
              data = hawks, Hess=TRUE)


reg.2x <- na.omit(hawks[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric +
                   married_dummy + province, 
                 data=hawks)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 26:  Insecurity and Support for Female Leadership (Doves)

doves <- data2[data2$Taliban_peace_support == "Strongly support", ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province,
            data = doves) 


reg.1x <- na.omit(doves[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
                female + age + edu_numeric +
                married_dummy + province,
              data = doves, Hess=TRUE)


reg.2x <- na.omit(doves[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric +
                   married_dummy + province, 
                 data=doves)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



##### 6. H2: Insecurity Decreases Support for Female Leadership (SI Section 3.4) ######

 

## Female Respondents Only ##


women <- data2[data2$female == 1, ]


## SI Table 27: Insecurity and Support for Female Leadership (Female Respondents Only)

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = women) 


reg.1x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
                + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = women, Hess=TRUE)


reg.2x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=women)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 28: Fear and Support for Female Leadership (Female Respondents Only)


reg.1 <- lm(female_support_continuous ~ fear_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = women) 


reg.1x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


#  ordinal

reg.2 <- polr(female_support_ordered  ~ fear_dummy + 
               + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = women, Hess=TRUE)


reg.2x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ fear_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=women)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



## SI Table 29: Insecurity and Support for Female Leadership 
## (Female Respondents, Remove wave fixed effects)


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province,
            data = women) 


reg.1x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = women, Hess=TRUE)


reg.2x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=women)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 30: Insecurity and Support for Female Leadership
## (Female Respondents, Balkh, Kunduz and Sar-e-Pul)


women2 <- data3[data3$female == 1, ]



reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy  + province + wave_factor,
            data = women2) 


reg.1x <- na.omit(women2[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal - province can't run with this.

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               age + edu_numeric + 
               married_dummy + wave_factor,
             data = women2, Hess=TRUE)


reg.2x <- na.omit(women2[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=women2)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 31: Insecurity and Support for Female Leadership 
## (Female Respondents, Hawks)

women_hawk <- women[women$Taliban_peace_support != "Strongly support", ]

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province ,
            data = women_hawk) 


reg.1x <- na.omit(women_hawk[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = women_hawk, Hess=TRUE)


reg.2x <- na.omit(women_hawk[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=women_hawk)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 32: Insecurity and Support for Female Leadership
## (Female Respondents, Doves)

women_dove <- women[women$Taliban_peace_support == "Strongly support", ]

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province,
            data = women_dove) 


reg.1x <- na.omit(women_dove[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = women_dove, Hess=TRUE)


reg.2x <- na.omit(women_dove[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=women_dove)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## H2: Insecurity and Support for Female Leadership Among Male Respondents
## (SI Section 3.4.2)



## SI Table 33: Insecurity and Support for Female Leadership (H2)

men <- data2[data2$female == 0, ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = men) 


reg.1x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = men, Hess=TRUE)


reg.2x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=men)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



## SI Table 34: Fear and Support for Female Leadership (H2)

reg.1 <- lm(female_support_continuous ~ fear_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = men) 


reg.1x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ fear_dummy + 
               + age + edu_numeric +
               married_dummy + province + wave_factor,
             data = men, Hess=TRUE)


reg.2x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ fear_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=men)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 35: Insecurity and Support for Female Leadership
## (Male Respondents, No Wave Fixed Effects)


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province,
            data = men) 


reg.1x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


#  ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = men, Hess=TRUE)


reg.2x <- na.omit(men[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=men)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 36: Insecurity and Support for Female Leadership
## (Male Respondents; Balkh, Kunduz and Sar-e-pul)


men2 <- data3[data3$female == 0, ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = men2) 


reg.1x <- na.omit(men2[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               age + edu_numeric +
               married_dummy  + wave_factor,
             data = men2, Hess=TRUE)


reg.2x <- na.omit(men2[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=men2)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 37: Insecurity and Support for Female Leadership 
## (Male Respondents, Hawks)


men_hawk <- men[men$Taliban_peace_support != "Strongly support",]

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province,
            data = men_hawk) 


reg.1x <- na.omit(men_hawk[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1


#  ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = men_hawk, Hess=TRUE)


reg.2x <- na.omit(men_hawk[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=men_hawk)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 38: Insecurity and Support for Female Leadership 
## (Male Respondents, Doves)

men_dove <- men[men$Taliban_peace_support == "Strongly support",]

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province,
            data = men_dove) 


reg.1x <- na.omit(men_dove[ , c("province", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$province))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               + age + edu_numeric +
               married_dummy + province,
             data = men_dove, Hess=TRUE)


reg.2x <- na.omit(men_dove[ , c("province", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$province))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province, 
                 data=men_dove)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## Results when subsetting to provinces where we did
## not run our survey experiment


## SI Table 39 Compare these results to non-experiment provinces

  
data3b <- data2[data2$experiment == 0,]


# H1

reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + female + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = data3b) 


reg.1x <- na.omit(data3b[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# See with ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               female + age + edu_numeric +
               married_dummy + wave_factor,
             data = data3b, Hess=TRUE)


reg.2x <- na.omit(data3b[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   female + age + edu_numeric + 
                   married_dummy + province + wave_factor, 
                 data=data3b)



stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



# H2 (Female Respondents)

## SI Table 40: Insecurity and Support for Female Leadership
## (Female Respondents, Non-Experimental Survey Districts)


women2b <- data3b[data3b$female == 1, ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy  + province + wave_factor,
            data = women2b) 


reg.1x <- na.omit(women2b[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1




# See with ordinal 

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
               age + edu_numeric + 
               married_dummy + wave_factor,
             data = women2b, Hess=TRUE)


reg.2x <- na.omit(women2b[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                   age + edu_numeric +
                   married_dummy + province + wave_factor, 
                 data=women2b)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)


## SI Table 41: Insecurity and Support for Female Leadership
## (Male Respondents, Non-Experiment Provinces)

men2b <- data3b[data3b$female == 0, ]


reg.1 <- lm(female_support_continuous ~ violence_experience_dummy + age + edu_numeric +
              married_dummy + province + wave_factor,
            data = men2b) 


reg.1x <- na.omit(men2b[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1


# ordinal

reg.2 <- polr(female_support_ordered  ~ violence_experience_dummy + 
                age + edu_numeric +
                married_dummy  + wave_factor,
              data = men2b, Hess=TRUE)


reg.2x <- na.omit(men2b[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


# use multinomial logit (base is equal)


reg.3 <- multinom(female_support_unordered ~ violence_experience_dummy +
                    age + edu_numeric +
                    married_dummy + province + wave_factor, 
                  data=men2b)


stargazer(clreg1, clreg2, reg.3)
stargazer(reg.1, reg.2, reg.3)



##### 7: Mechanisms #######


## SI Table 42

reg.1 <- lm(prov_gov_support_rank ~ violence_experience_dummy + age + edu_numeric +
               married_dummy  + province + wave_factor,
             data = women)

reg.1x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 

clreg1 # negative relationship


reg.2 <- glm(prov_gov_support_dummy ~ violence_experience_dummy + age + edu_numeric +
              married_dummy  + province + wave_factor,
            family=binomial,
            data = women)

reg.2x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 

# negative relationship


stargazer(clreg1, clreg2)
stargazer(reg.1, reg.2)



## SI Table 43

# Are you opposed to women representing you in prov gov? 1 = YES -> women less supportive
# This question was only asked in earlier rounds

reg.1 <- lm(women_prov_binary ~ violence_experience_dummy + age + edu_numeric +
               married_dummy  + province + wave_factor,
             data = women)

reg.1x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.1)))]) 

clreg1 <- coeftest(reg.1, vcov=vcovCL(reg.1, factor(reg.1x$wave_factor))) 


reg.2 <- glm(women_prov_binary ~ violence_experience_dummy + age + edu_numeric +
               married_dummy  + province + wave_factor,
             family=binomial,
             data = women)

reg.2x <- na.omit(women[ , c("wave_factor", all.vars(formula(reg.2)))]) 

clreg2 <- coeftest(reg.2, vcov=vcovCL(reg.2, factor(reg.2x$wave_factor))) 


stargazer(clreg1, clreg2)
stargazer(reg.1, reg.2)


